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We report the first measurements of phonon dispersion curves on the (001) 
surface of the strong three-dimensional topological insulator Bi2Se3. The sur- 
face phonon measurements were carried out with the aid of coherent helium 
beam surface scattering techniques. The results reveal a prominent signature 
of the exotic metallic Dirac fermion quasi-particles, including a strong Kohn 
anomaly. The signature is manifest in a low energy isotropic convex dispersive 
surface phonon branch with a frequency maximum of 1.8 THz, and having 
a V-shaped minimum at approximately 2kF that defines the Kohn anomaly. 
Theoretical analysis attributes this dispersive profile to the renormalization 
of the surface phonon excitations by the surface Dirac fermions. The contri- 
bution of the Dirac fermions to this renormalization is derived in terms of a 
Coulomb-type perturbation model. 
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The recent discovery of the new class of materials coined topological insulators (TIs) [HI 
El El m m has attracted much interest and excitement. Strong TIs have exotic metallic surface 
states protected by time-reversal invariance (TRI). The presence of the metallic boundary is 
dictated by the recent Z2 classification of TRI insulators into an ordinary insulator class (Z2- 
even), including the vacuum, and a topological insulator class (Z2-odd). Members of each class 
can be adiabatically converted to each other, but not into members of the other class [l6l|71[8l. 
Ordinary and topological phases are separated by a topological phase transition, where the bulk 
gap has to vanish at the transition point. Thus, the ensuing scenario for a strong TI depicts 
a topological bulk (Z2-odd), embedded in vacuum (Z2-even), a case which demands that the 
boundary (surface) be gapless. This gives rise to surface states involving massless Dirac fermion 
quasi-particles (DFQs) having an odd number of Dirac cones in the surface Brillouin zone 
(SBZ). Moreover, the strong spin-orbit coupling leads to a definite helicity whereby the spin 
is locked normal to the wavevector of the electronic state. A fundamental constraining feature 
of such spin-textured surface states is their robustness against spin-independent scattering, an 
attribute that protects them from backscattering and localization [|9l[T0l. 

In this paper we address, experimentally and theoretically, a manifestation of the exotic 
surface DFQs that has received little attention: their response to surface phonon excitations, 
namely surface electron-phonon scattering and the ensuing phonon energy renormalization. In 
this work, we employ elastic and inelastic helium atom surface scattering (HASS) techniques to 
determine the surface structure and obtain the surface phonon dispersion curves, respectively. 
It is well established that He atoms are scattered by the oscillations of the surface electron 
density about 2-3 A away from the first atomic layer, and thus HASS is very sensitive to 
phonon-induced surface charge density (SCD) oscillations, even those induced by subsurface 
second-layer displacements [fm[T2l . Consequently, HASS intensities carry direct information 
on the SCD oscillations associated with surface phonons and ultimately on the surface electron- 
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phonon (e-p) interaction. HASS, therefore, is an ideal technique to investigate the collective 
Dirac electron states response to surface lattice excitations. 

Because the family of strong TIs 612X3 (X=Se, Te) was found to have a single Dirac cone 
with a Dirac point at the SBZ center [fT3l [T4l [T5l [T6l [TTl . it has been the focus of ongoing 
intense experimental and theoretical studies. Angle- and spin-resolved photoemission spec- 
troscopy measurements provided strong evidence of the existence of the linearly dispersive 
Dirac cone and the spin texture. Moreover, scanning-tunneling microscopy confirmed the ab- 
sence of backscattering, namely, scattering between states of opposite momentum and opposite 
spin [HI [181. We selected to study Bi2Se3 because of its simplicity, its relatively wide gap of 300 

meV, and the rich trove of information available about its DFQ surface states lfT9l[20l[2n[22l[23ll . 

• • • 




• • • 

Figure 1: The extended surface Brillouin zone, showing high-symmetry points F, M and K. The F-M 
direction lies along q^;, while the F-K-M line is traced along qj,. The GjS are surface reciprocal lattice 
vectors. 



The (001) Bi2Se3 crystal surface belongs to the pQmm two-dimensional space group. The 
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corresponding SBZ is shown in figure [l| with high-symmetry points F, M and K. Extensive 
elastic and inelastic He scattering measurements of the Bi2Se3(001) surface were carried out 
[|24ll . The processed inelastic data is presented in figure [s] along the F-M, F-K and K-M direc- 
tions as orange dots together with their respective error bars. 

In order to interpret and fit the experimental inelastic data, we carried out phenomenological 
surface lattice dynamical calculations, based on the pseudo-charge model (PCM) [|25l l26l l24l 
and applied to slab geometries containing 30 quintuples. The physical role of the surface DFQ 
states in the phonon energy renormalization was subsequently studied with the aid of a micro- 
scopic model based on Coulomb-type perturbations. The model was incorporated in calcula- 
tions of density-density correlations that take into account the helicity and linear dispersion of 
the surface DFQ states. 

612863(001) surface phonon dispersions: Experimental and PCM calculations results 

The slab calculations were preceded by a bulk calculation [|24ll based on PCM and fitted to 
available Raman and IR data [|27ll28l . To obtain a best-fit to the measured surface phonon dis- 
persion curves the following changes to bulk parameter values were made: The surface Se-Bi 
force constant was reduced by 25% from its bulk value; this is a reasonable change since the 
non-metallic bonding in the two topmost layers is effectively reduced to allow for the emer- 
gence of the metallic electrons. A new planar force constant involving surface Se-Se ions was 
introduced. To account for the metallic deformability of the surface DFQs we reduced the mag- 
nitude of the pseudocharge-ion coupling constant Ts in the topmost pyramid involving Se and 
Bi from its bulk value Tl by about ^ = ~ 13% Physically, AT5 accounts for the 

^ B B 

extra screening provided by the DFQ surface states, which is proportional to the corresponding 
Fermi surface density of states, Vp oc Ep oc kp. Hence AT5 oc kp. 

However, to underscore the significance of the experimental results and their interpretation 
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in terms of PCM fitting, we start by presenting in figure [2] the surface plionon dispersions for 
a model employing the bulk value of the ion-pseudo-charge coupling parameter at the surface, 
namely Ts = = 8.07 N/m [|24ll . The calculated surface phonon dispersion curves are 
presented in figure [2] as black dots, and the gray background represents the projection of the 
bulk bands onto the SBZ. The important feature to be noted is the presence of a surface Rayleigh 
branch extending from cj = 0tO(X'^3.7 THz. 
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Figure 2: Surface phonon dispersion curves (black dots) along the F-M and F-K directions, for 

Ts = = 8.07 N/m bulk value. The gray background represents the projection of the bulk bands on 

the SBZ. 

By contrast, the calculated surface phonon dispersion curves that provide a best-fit for the 
experimental data, correspond to a reduced value of Ts of 7.05 N/m. The calculated surface 
phonon dispersion curves are shown in figure[3]as black dots, superimposed on the experimental 
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Figure 3: Surface phonon dispersion curves along the F-M and F-K directions: Experimental data 
appear as orange dots with error bars reflecting instrument resolution, while the calculated dispersions 
curves, using PCM with Ts = 7.05 N/m, are represented by black dots. The gray background represents 
the projection of the bulk bands on the SBZ. The red dot on F at 160 cm^^ represents an experimental 
surface Raman mode reported in Ref. ll29ll . 

data which is displayed as solid orange dots with error bars. The single red dot at the F-point, 
at about 5 THz, is a Raman surface frequency reported in ||29l . 

We shall focus our attention on two important observations regarding figure |3j First, we 
should discern the absence of a surface Rayleigh branch in both experimental and calculated 
dispersion curves. Second, we observe the emergence of a "prominent" isotropic parabolic 
dispersion branch centered at the F-point with a frequency of 1.8 THz (7.4 meV); the PCM 
calculations confirm that at F it has a vertical shear (z) polarization, with M|e/^Bi = 2.8, 
where is the respective ionic displacement. This branch terminates in a V-shaped minimum 
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located at g ~ 0.2 A^^, a value that roughly corresponds to 2kp of the DFQs, and thus signals 
the manifestation of a strong Kohn anomaly [|30l . The isotropy of this branch and its apparent 
termination at 2kp can be explained by a scenario involving the DFQ surface states, in particular 
their isotropic Fermi surface. In this scenario, the V-shaped feature marks the boundary between 
an operative DFQ screening for q < 2kF, and its suppression above this value, which is a typical 
signature of a Kohn anomaly. Lattice dynamics calculations reveal some bulk penetration of 
vertical shear modes for q > 2kp reflecting a diminished role of DFQ screening and more 
compatibility with the insulating bulk. We shall establish below the intimate hnk between the 
dispersive character of this branch and the surface DFQ states response to ionic displacements. 

The polarization and other properties of the remaining surface phonon dispersive branches 
are discussed in [|24l . 

Perturbative Coulomb interactions results 

We now briefly describe the microscopic model used to study the contributions of the surface 
Dirac electronic states that define the dispersive character of the prominent phonon branch. Our 
experimental results demonstrate that this surface phonon branch is "optical" in nature, and thus 
should be derived from a Coulomb-type perturbative mechanism. Consequently, we construct 
a perturbative Hamiltonian describing the interaction of the surface Dirac electrons with the 
electric field resulting from ionic displacements, namely, a linear coupling of the lattice ionic 
displacement to the DFQ density of the form 



where Pci(i') is the 2-dimensional electron density, and is the displacement of the j^^ surface 
ion whose equilibrium position is 'R.f^- A is a vector accounting for the position-dependent cou- 
pling of the charge density to ionic displacements. As shown in [|24ll . after second quantization. 




(1) 
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with 6q ,y and Ck, a being the annihilation operators for a surface phonon mode (q, 7) with po- 
larization index 7 and an electron with momentum k and spin a, respectively, the Hamiltonian 
takes the form 



"Hcl-ph — gq,7 Ck_|_q^ „ Ck, a ^q,7 , (2) 

where Aq^^ = (6q,^ + &q,^)- The electron-phonon coupling constant is 



Nh ^ ^ . Nh ^ 

A is the surface area, the number of primitive cells in A, M the ionic mass. O0q)y is the bare 
phonon frequency of mode (q, 7) and e^(q) the corresponding polarization vector. Aq denotes 
the Fourier transform of A(r). 

A detailed derivation given in [|24ll leads to the Dyson equation 

where cuq^j is the renormalized surface phonon frequency, and 21 is the surface primitive cell 
area. 11 and e are the polarization and dielectric functions in the random phase approximation, 
respectively, expressed in the helicity basis; the former can be decomposed into two contribu- 
tions: one due to intra-band and the other due to inter-band excitations as shown in figure |4j 



n(q,u;) = n-'^^(q,a;) + n'"'-(q,u;). (5) 
Explicit expressions for H'""^^ and H'"''^'^ are given in [|24]| . 



In order to solve (43 1 self consistently, we need to specify the momentum dependence of the 
coupling function Aq,^. We argue that in the range of |q| < 2k p, the ionic screened potential 
V{q) has very weak dependence on q, since 2kp < q^F = ilj^pp — O-^ '^he Thomas- 
Fermi wave vector. Moreover, invoking the sagittal-plane symmetry classification of the surface 
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Figure 4: Intra- and Inter-band transitions of DFQs that contribute to the renormalization of the promi- 
nent surface phonon branch, q is the phonon wave vector. 



phonon modes, we specify the general polarization vector of the even parity mode as 

e^(q) ^ e_L(q) ~ z, e||(q)~q. (6) 

ex(q) accounts for the vertical shear component (ionic oscillations in a direction normal to the 
surface plane) while e|| (q) is the longitudinal component. The corresponding coupling constant 
then reads 

Aq,7 = Aq • e^(q) = Aq-e^(q) + Aq-e|j(q) 

= A^(q) + A||(q) . (7) 

In view of the near constancy of V^(q) for q < 2kp, and the fact that the electron-phonon 
coupling involves the gradient of the screened potential, we write 

Aq = Af + ^a[°) (8) 

where go will be appropriately chosen as 2kF. Notice that at F (q = 0) the polarization is pure 
vertical sheer, in agreement with the surface symmetry and the PCM. 
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Figure 5: Renormalized surface topological phonon dispersion curve, superimposed on the correspond- 
ing experimental data. 



The Dyson equation now becomes 



to first order in q. 

The model depends on two parameters: the bare phonon frequency and the Fermi wave 
vector. We identified the former with the experimental value of (X'(q = 0) = 1.8 THz, where 
the DFQ response vanishes, kp was derived from a sample carrier concentration of —1.9 x 
10^^/cm'^, obtained from Hall measurements, which correspond to a Fermi energy of about 
300 meV and a Fermi wavevector of kp = 0.1 A^^, which is consistent with previous values 
reported fromphotoemission measurements [[T9ll2Tll22l . The coupling parameters A|^^ and A||°^ 
were left as fitting variables. Solutions for a;(q) were obtained by carrying out the integrals for 



n(q, w), and solving Eq. (43) iteratively. The calculated best-fit dispersion curve is shown in 
figure [5| superimposed on the experimental data. The agreement is quite good. 
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The best fit parameters are 




{X^y = 10^ (meV) 




A|| 



= 0.65 . 



(10) 



With 



fi2 

MS21 



4 X 10 ^ meV, we obtain 



= 50 ey • A . 



This yields a real space value of 3.4 eV/A, which is quite reasonable since it falls in the range 



of typical Coulombic interactions. 

Finally, we demonstrated above that a TI actually presents a composite system consisting 
of an "ultra-thin metallic film" and an underlying insulating substrate. As such, the phonons 
in the metallic film can hardly penetrate into the substrate bulk in order to establish an acous- 
tic Rayleigh branch. However, the unique feature of this system is that for g > 2kF, as DFQ 
screening becomes gradually suppressed, the surface film becomes almost homogeneous with 
the underlying substrate, and the corresponding modes gradually penetrate into the bulk, ush- 
ering the V-shaped dispersion. We designate this unique behavior as a strong Kohn anomaly. 
Moreover, we should emphasize here that theoretical modeling and analysis of the interaction 
between the surface phonons and the DFQs demonstrate that the linear dispersion and isotropy 
of the DFQs are responsible for the profile of the prominent surface phonon branch for q < 2kF- 

Methods 

Theoretical Calculations 

The pseudo-charge model (PCM) was incorporated into a slab-geometry program to calculate 
the surface phonon dispersion curves. PCM includes the electronic degrees of freedom as mass- 
less particles (adiabatic approximation) with dynamical variables associated with symmetry- 
adapted deformation functions. Simple eUmination of these dynamical variables leads to an 
effective screening of the ionic motions. 
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Experimental Methods 

Bi2Se3 single crystals were prepared using the vertical Bridgman method. The initial stoichio- 
metric mixtures of Bi and Se powder of 5N purity were sealed in an evacuated conical quartz 
ampoule. Typical crystal samples used were about 3mm x 3mm x 2mm in size, and were cleaved 
in situ under ultra-high vacuum conditions, with the aid of an attached post. Inelastic HASS 
measurements were carried out with an angle-resolved detector equipped with a time-of-flight 
(TOF) system using a gating technique based on electronic excitations of helium atoms to the 
2^S triplet metastable state. Data processing of TOF spectra was done with the aid of scan 
curves, namely AE{AK), where AE and AK are the energy exchange and momentum trans- 
fer, respectively. 

A more detailed account of the Methods is given in the supplementary material. 
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Supplementary Material 

1. The topological insulator Bi2Se3 




Figure 6: Crystal Structure of Bi2Se3. (Left) A unit cell showing the quintuple structure. (Right) A top 
view showing the stacking arrangement. 



Bi2Se3 has a rhombohedral crystal structure with space group Dj^ (R3m). The primitive 
cell, shown in figure|6} contains a binary axis (with twofold rotation symmetry), a bisectrix axis 
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Figure 7: Brillouin Zones (BZ). (a) A schematic view of the bulk three-dimensional BZ of Bi2Se3 and 
the two-dimensional BZ of the projected (001) surface, (b) The extended surface BZ with reciprocal 
lattice vectors Gi and G2. Experimental measurements were carried out along the two high symmetry 
directions, F-M and F-K, which are aligned along qx and qy, respectively. 



(appearing in the reflection plane) and a trigonal axis (with threefold rotation symmetry). The 
primitive translation vectors in the rhombohedral and hexagonal basis are 



~ [ 2' 6 ' 3)' ^ ~ [2' 6 ' ' 3 - (^^' 6 ' ' 
tf = (a, 0, 0),t^ = (^|, 0^ , = (0, 0, c) , (11) 

where a and c are lattice constants of the hexagonal cell. The corresponding reciprocal vectors 
are given as 

= (-,, -4, b? = fi, -4, ^\ b? = fo, 



c I a \ b c I a \ b c j a 

b« = ( 1. 0) -, b» = (a. ^, 0) -, b» = (0. 0, 1) (12) 
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The bulk Brillouin zone and surface Brillouin zone (SBZ) are shown in figure |7ja; while an 
extended version of the SBZ is displayed in figure |7jb. Bi2Se3 displays a layered structure 
consisting of five atomic layers as a basic unit (cell), ordered in the Se(l)-Bi-Se(2)-Bi-Se(l) 
sequence, designated a quintuple layer (QL), as shown in figure[6} Se(2) is a center of inversion. 
Within each QL, the inter-layer bonding is strong and exhibits a predominant covalent character. 
The inter-quintuple bonding is of the van der Waals-type. 

2. Experimental Setup and Procedures 
A. Crystal Preparation 

Bi2Se3 single crystals were prepared using the vertical Bridgman method. The initial stoichio- 
metric mixtures of Bi and Se powder of 5N purity were sealed in a conical quartz ampoule and 
evacuated to 10 mTorr after repeated Argon gas purging. Preliminary synthesis and homoge- 
nization was carried out in a horizontal tube furnace at 650 C for 48 hours. The sealed ampoules 
were then passed through a vertical Bridgman furnace with a hot zone at 850 C and a thermal 
gradient of 2 C/cm near the solidification point 705 C. The pulling rate was kept at 0.2 mm/hr 
after soaking at 850 C for 24 hours. The grown single crystals were 3 cm long and about 1.2 
cm in diameter, with good optical quality. They are easy to cleave with crystal planes perpen- 
dicular to the trigonal c-axis. Hall conductivity measurements on these samples gave carrier 
concentration of about — 1.9 x 

Typical crystal samples used were about 3mm x 3mm x 2mm in size, with its exposed surface 
parallel to the QL planes. The crystals were attached to an OFHC copper sample-holder by 
conductive silver epoxy. A cleaving (peeling) post was attached to the top sample surface in 
a similar way. The prepared sample holder was mounted on a sample manipulator equipped 
with XYZ motions as well as polar and azimuthal rotations. The pressure in the Ultra-High 
Vacuum (UHV) chamber was maintained at 10^° torr throughout the experiment to ensure the 



15 



cleanliness of the sample surface during measurement. In situ cleaving under UHV conditions 
was effected by knocking off the cleaving post. All measurements were performed with the 
sample surface at room temperature. Immediately after cleaving, the quality of the long-range 
ordering on the surface was confirmed by the appearance of sharp diffraction LEED spots. 

B. Helium beam monochromator 

A supersonic mono-energetic coUimated helium beam, with velocity resolution better than 
1.4 %, was generated by a nozzle- skimmer assembly and 2mm diameter coUimating slits. The 
average beam velocity was varied by attaching the nozzle reservoir to a closed-cycle helium 
refrigerator, and controlling the reservoir temperature, in the range 300 K to 110 K, with the 
aid of a digital temperature controller (Scientific Instruments Model 9700) and a diode sensor 
attached to the reservoir. This allows the beam energy to vary in the range 65 meV to 21 meV. 
Polar rotation of the sample was used to vary the incident angle 9i with respect to the surface 
normal, while the azimuthal rotation was employed to align the scattering plane along a high- 
symmetry surface crystallographic direction. 

C. Scattered beam detection technique 

The scattered He beam was collected by an angle-resolved detector mounted on a two-axis 
goniometer, which allows the scattered angle 6*/ to be varied independently from 9i [|3T1l . and 
allows in- and out-of the scattering-plane measurements. As shown in figure [8| the detector 
is comprised of an electron gun and a multichannel plate (MCP) electron multiplier [|32ll . The 
electron gun generates a well-coUimated, monoenergetic electron beam crossing the He beam 
at right angles. The energy of the electron beam is tuned to excite the He atoms to their first 
excited metastable state 2^S (He*) upon impact. Deexcitation of a He* atom at the surface of 
the MCP leads to the ejection of an electron (Penning ionization) which generates an electron 
cascade that is collected by the anode of the multiplier. By electronically pulsing the electron 
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Figure 8: Schematic of the He beam scattering geometry, showing the metastable He* atom detection 
scheme. 

gun, a gate function is created for time-of-flight (TOF) measurements in the inelastic HASS 
mode. The details of the detection scheme are given in Ref [|32ll . 

Inelastic scattering measurements were carried out for in- and out-of- the- scattering-plane 
geometry. By writing the He-atom wave vector as k = (K, kz), where K is the component par- 
allel to the surface, conservation of momentum and energy for in-the-scattering-plane geometry 
can be expressed as 



AK = G + Q 

AE = hu}{Q) 



k f sin Of — ki sin 6i 



h^kj 



h^kf 



2MHe 2MHe 



(13) 
(14) 
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where subscripts i and / denote incident and scattered beams, respectively, and AK is the 
momentum transfer parallel to the surface. G is a surface reciprocal-lattice vector, Q is the 
surface phonon wave vector, hio^Q) is the corresponding surface phonon energy and Mhe is the 
mass of a He atom. By eliminating k/ from the above equations, one obtains the so-called scan 
curve relations which are the locus of all the allowed AK and AE as dictated by the geometry 
and the conservation relations. 



where Ei = h^kf /2Mue- The intersections of these scan curves with the phonon dispersion 
curves define the kinematically allowed inelastic events for a fixed geometric arrangement. 
Thus, by systematically changing E^, 6i, and 6f, the entire set of dispersion curves can be 
constructed. 

D. Typical Experimental Results 

Figure [o] (a) and (b) show typical diffraction patterns along the F-M and F-K directions, con- 
firming the hexagonal geometry of the surface and the lattice constant a = 4.15 A. Figure |9](c) 
and (d) show typical TOF spectra displaying diffusive elastic as well as inelastic peaks; Gaus- 
sian fits to the displayed peaks are shown in blue for the diffuse elastic peak, and in green for 
the inelastic ones. Figure |9](e) and (f) display typical scan curves with experimental inelastic 
events shown as red dots. 



In order to identify the character of the measured phonon dispersion curves in terms of irre- 
ducible representations (Irreps) and polarization vectors, we use an empirical lattice-dynamical 
model, the pseudo-charge model, that includes both Born- von Karman type direct ion-ion force 
constants, as well as indirect ion-ion adiabatic coupling through the mediating electrons. This 




(15) 



3. Pseudo-charge Phonon Model 11251. 1261. fl2l| 
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Figure 9: Typical experimental results from HASS measurements, a-b, Elastic diffraction patterns 
along F-M (a), and F-K (b) directions, c-d, Typical energy loss (gain) spectra from inelastic TOP 
measurements. Experimental data (circles) are resolved into inelastic peaks (green) and diffuse elastic 
peaks (blue), e-f, Scan curves with TOE experimental results (red dots). Each red dot represents a surface 
phonon event, and also appears as a green peak in c-d. 
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model has been particularly successful in reproducing both the surface dispersion curves and 
the HASS scattering amplitudes, derived from calculated surface charge deformations. 

The model accounts explicitly for the electronic degrees of freedom in the following man- 
ner: The electron density within each primitive cell i is expanded in terms of symmetry- 
adapted (SA) multipole components around selected Wyckoff symmetry points R^j, so that the 
corresponding pseudo-charge is written as 

ne{r) = ^ crfc lr,fc(r - R^^) (16) 

where F denotes an Irrep of the Wyckoff symmetry point-group and k indexes its rows; Ypk is a 
SA harmonic function. The expansion coefficients crk are treated as bona-fide time-dependent 
dynamical variables, alongside the ionic displacements that they couple to, namely 

crk{t) = + Acrfc(t) 

The coupling is treated within the adiabatic approximation. 

A Taylor expansion of the total potential energy in terms of ionic displacements Mq, 
where a and k denote Cartesian directions and sublattices, respectively, and pseudo-charge 
fluctuations Acrfe is given in matrix form by 

(17) 



V{r) = Vo + l- u u + (u • T . Ac + h.c.) + Ac • H • Ac 



2 

where u and Ac are vector matrices representing ionic displacements and charge deformations, 

respectively. $ is the ion-ion force-constant matrix, T the electron-ion coupling matrix and H 

is the electron-electron coupling matrix. The elements of these matrices are treated as empirical 

parameters to be fitted to experiment. 

The corresponding equations of motion are 

M,M£k) = - ^ (l I) up {^'k') - T. ([ Acrfc (/j) (18) 
mvAcrkm = - V T« f,) Ua {£' n') - Hr,km Acrk{ij) (19) 

a V / rv / 
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Figure 10: A schematic diagram of the lattice structure and the pseudo-charge model. Two quintuples 
(Sel-Bi-Se2-Bi-Sel) are shown in the diagram. The figure indicates that the surface pseudo-charges are 
more spread and easier to deform than their bulk counterparts. The figures to the right, show the dipolar 
pseudo-charge deformations for lattice distortions along the z direction. 



where and nir are the mass of ion and the electronic effective mass, respectively. 
Invoking the adiabatic approximation, mr = 0, we obtain 

Ac = H^T^u 

while translation invariance against an arbitrary rigid displacement Uq leads to 

Ac = H-^T^Uo 



and 







$ + TH^T^ 



Uo 



so that the effective ion self force constant matrix is expressed as 



K K 



I'k' ^ ' j,Tk ^ 



H 



-1 T-iT 

• -"-rfc 



j K 



(20) 



(21) 



(22) 



(23) 
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The Fourier transformed equation of motion reduces to 

a;^(q)u(q) = X>(q) u(q) 

2?(q) = ($(q) - T(q)H-i(q)Tt(q)) (24) 

where X'(q) is the dynamical matrix. 

We found that the 6c Wyckoff positions of the R3m space group were the most appro- 
priate to use as centers of pseudo-charge SA multipole expansion. These Wyckoff positions 
are identified as having coordinates (0, 0, ±z) and Ci,^ point-group symmetry; they define the 



axes of the tetrahedral pyramids shown in figure 10 Among these points, the center of ex 



pansion was chosen to be the pyramid center, again, depicted in figure 10 has h^reps 
Ai (with SA dipolar harmonic z) and E (with SA dipolar harmonics x, y). In order to minimize 
the number of empirical constants employed, we opted to include only the Ai SA fluctuations. 



The right part of figure 10 shows a schematic representation of the pseudo-charge deformation 



associated with Ai symmetry. 

A. Bulk phonon dispersion curves 

The initial fitting procedure was carried out for bulk phonons where Raman and IR phonon 
frequencies are available in the literature [|27ll28l . To our knowledge no neutron scattering data 
has been reported for Bi2Se3. Central ion-ion interaction potentials, v{r), were adopted; with 
force-constant matrix elements given by 

Xq, / Xq, X 1 



where 



A 

Q^2 



, and B = ^ 

r=ro or 
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r=ro 
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Figure 1 1 : Calculated phonon dispersion curves for the bulk of Bi2Se3 along the T-Z direction, based 
on the pseudo-charge model. The experimental data are from Ref.|27 |. The red dots correspond to exper- 
imental Raman frequencies, and the blue circles correspond to experimental infrared (IR) frequencies. 



where ro is the equilibrium bond length. A third parameter, T, defines the ion-pseudocharge 
coupling. We used a value of Ti for all pyramids involving Sel-Bi and T2 for those involving 
Se2-Bi. In the present model, we neglect pseudocharge-pseudocharge interactions, thus render- 
ing the matrix H diagonal; its elements are determined from translation invariance constraints. 



The best-fit bulk phonon dispersion curves are presented in figure 1 1 along the high-symmetry 
direction T-Z. Raman and IR experimental data are depicted as solid red dots and open cir- 
cles, respectively. The corresponding fitting parameters are given in Table 1, where the ion-ion 
parameters are specified for each bond-type; the Sel-Sel bond being the van der Waals inter- 
quintuple one. We also present both experimental frequencies and their calculated counterparts 
in Table 1 for comparison. For reference, the bulk BZ is shown in figure |7]-a. 
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Table 1: Best-fit parameters (in N/m) for bulk pseudocharge model, together with experimental 
(Ref.[|27ll28|) and calculated frequencies. 



Ion-Ion Interaction 


lon-Pseudocharge Interaction 


Frequencies at F-point (THz) 


Bond A B 


Position Constant 


Calculated Experiment 


Sel-Sel 1.0 0.1 
Sel-Bi 13.5 1.35 
Bi-Se2 3.0 0.3 
Bi-Bi 2.0 0.2 


(Bi-Sel) 8.07 
T| (Bi-Se2) 7.46 


Ai„ 5.496 NA 
Aig 5.235 5.235 
Aij, 3.554 NA 
Aig 2.160 2.16 
Eg 3.979 3.945 
E„ 3.945 3.87-4.02 
E„ 1.836 1.83-1.95 
Eg 1.094 NA 



B. Surface phonon dispersion curves 

Lattice-dynamics calculations of the surface phonon dispersions based on a pseudocharge lattice- 
dynamical model [|26l IT2l and a slab geometry containing 30 quintuples were carried out. The 
(001) Bi2Se3 crystal surface belongs to the p6mm two-dimensional space group. The corre- 
sponding surface Brillouin zone (SBZ) is shown in figure |7} The high-symmetry points are 
r, M and K. The F-M direction lies along q^., while the F-K-M line is traced along q^. 

To obtain a best-fit to the measured surface phonon dispersion curves the following model 
changes to bulk parameter values were made: The surface Se-Bi force constant was reduced by 
25% from its bulk value; this is a reasonable change since the non-metallic bonding in the two 
topmost layers is effectively reduced to allow for the emergence of the metallic electrons. Also 
a new planar force constant involving surface Se-Se ions was introduced. To account for the 
metallic deformability of the surface DFQs we reduced the magnitude of the pseudocharge-ion 
coupling constant in the topmost pyramids involving Se and Bi by about = '^^r^i^^ — 13%. 
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Physically, ATs accounts for the extra screening provided by the DFQ surface states which 
is proportional to the corresponding Fermi surface density of states Vp oc Ep oc kp. Hence 
ATs oc kp. Table |2] lists all the modified surface parameters. 

However, to underscore the significance of the experimental results and their interpretation 
in terms of model fitting, we start by presenting the surface phonon dispersions for a model 
employing the bulk value of the ion-pseudo-charge coupling parameter at the surface, 
figure 12 The calculated surface phonon dispersion curves are presented in figure [2] as black 



dots, and the gray background represents the projection of the bulk bands onto the SBZ. The 
important feature to be noted is the presence of a surface Rayleigh branch extending from a; = 

25 



to w ~ 3.7 THz. 



Table 2: Best-fit surface parameters (in N/m) of the pseudo-charge slab model calculation. 



Surface Ion-Ion Interaction Surface lon-Pseudocharge Interaction 



Ion-Ion 


A 


B 


Position 


Constant 


Sel-Sel 
Sel-Bi 


2.0 
10.0 


0.2 
1.0 


T5(Bi-Sel) 
ATs = T^-Ts=im 


7.05 



M K r M 




1.4 1.2 1.0 0.8 0.6 0.4 0.2 /-> 0.2 0.4 0.6 0.8 



Figure 13: Surface phonon dispersion curves along the F-M and F-K directions: Experimental data 
appear as orange dots with error bars reflecting instrument resolution, while the calculated dispersions 
curves, using PCM with Ts = 7.05 N/m, are represented by black dots. The gray background represents 
the projection of the bulk bands on the SBZ. The red dot on F at 160 cm~^ represents an experimental 
surface Raman mode reported in Ref. li29il . 



By contrast, the calculated surface phonon dispersion curves that provide a best-fit for the 
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experimental data, correspond to a reduced value of r5 of 7.05 N/m. The calculated best-fit 



surface phonon dispersion curves are presented in figure 13 as black dots, superimposed on the 
experimental data which is displayed as solid orange dots with error bars. The single red dot at 
the F-point, at about 5 THz, is a Raman surface frequency reported in [|29ll . 

We focus our attention on two important observations regarding figure [3} First, we should 
discern the absence of a surface Rayleigh branch in both experimental and calculated dispersion 
curves. Second, we observe the emergence of a "prominent" dispersion branch centered at the 
F-point with a frequency of 1.8 THz (7.4 meV); the model calculations confirm that it has a 
vertical shear (z) polarization at the F-point, with u^sel'^^Bi — 2-8. This branch exhibits an 
anomalous isotropic parabolic dispersion, centered at F, and terminates in a V-shaped minimum 
located at g ~ 0.2 A, a value that roughly corresponds to 2kF of the DFQs, and thus signals 
the manifestation of a strong Kohn anomaly [|30l . The isotropy of this branch and its apparent 
termination at 2/c^ can be explained by a scenario involving the DFQ surface states, in particular 
their isotropic Fermi surface. In this scenario, the V-shaped feature marks the boundary between 
an operative DFQ screening for q < 2kp, and its suppression above this value, which is a typical 
signature of a Kohn anomaly. Lattice dynamics calculations reveal some bulk penetration of 
vertical shear modes for q > 2kF reflecting a diminished role of DFQ screening and more 
compatibility with the insulating bulk. We shall establish below the intimate link between the 
dispersive character of this branch and the surface DFQ states response to ionic displacements. 

4. Hamiltonian of the electron-phonon interacting system 

In this section we discuss the renormalization of the surface phonon's frequencies on the (001) 
surface of Bi2Se3 due to the electron (DFQ)-phonon interaction in terms of a microscopic model 
that incorporates the surface Dirac Hamiltonian, the non-interacting surface phonon Hamilto- 
nian and the DFQ-surface phonon interaction Hamiltonian. The latter involves linear coupling 
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of the lattice ionic displacement to the DFQ density. As we shall demonstrate below, the model 
provides a very good fit to the experimental results for the "prominent" dispersion branch that 
has been Unked to DFQ screening above. The model also establishes the nature of the screening 
mechanism to surface phonon excitations, and, thus, complements the analysis provided by the 
pseudo-charge model. 

The system is described by the Hamiltonian 

"H = "Hel + ^ph + "Hel-ph = "Ho + "Hel-ph , (25) 

whereby Ho comprises the non-interacting part of the Hamiltonian while Hei-ph describes the 
interacting part, all of which will be specified hereafter. 

We take the (001) surface of Bi2Se3 to be the x-y plane, with area A normal to the unit 
vector z, and consider a two-dimensional lattice containing N ions, each of mass M, which 
are labeled by the index j, and assign to each of these atoms an equilibrium position Rj°^ A 
displacement of the j-th ion about its equilibrium position R^°^ can be expanded as 

q,7 

where Iq^^ = j^^m e^(q) is the polarization vector. The harmonic modes of vibration 
are then described by the standard phonon (lattice) Hamiltonian 

q,7 ^ / 

where one identifies as the creation operator of a phonon with energy Q^lf = h ujq)y and 
polarization index 7. 

The electronic surface states of Bi2Se3 form a two dimensional Dirac metal, whose low 
energy physics is well described by the Hamiltonian 

^el = ^ 4 z • (kx 0-) - //] V'k , (28) 

k 
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where tpk = i'^^^ ] is the two-component electron spinor operator at wave vector k, Vi? is the 
Fermi velocity, fi is the Fermi energy (which lies above the Dirac point) and cr = (cti, 0-2) is 



the vector containing the first two Pauli matrices. The Dirac Hamiltonian (|28J) is diagonal in the 
helicity basis ^'k ~ ' 



7k. 

= t^k^k 

1 fie^^- 1\ ^ fky\ (29) 



yielding 



nei = J2Yl (7k)^ 7k, e£ = a hvp \k\-fi. (30) 

k a=± 

We consider an interaction between the electron density and the displacement (where 
has both in-plane and out-of-plane components) of the j"^ ion about its in-plane equilibrium 
position 'R.f^- The electron-phonon interaction can be generically written as 

TV 

(31) 



^ei-ph = / d^rpe^r) A(r - Rf ) • Mj . 

.7=1 



Here, Poi(i") = c)^(i')ccr(i') is the electron surface density operator and A(r — R^" ) is a 
position dependent vector function (with units of energy per length) characterizing the electron- 
phonon coupling. 



We now proceed to represent the interaction Hamiltonian (31) in Fourier space. This is 
achieved by expressing: 

q 

Pei(r) = 4(r)c.(r) = \ H ^ "''^''^ 4+c,.Ck,. . (33) 



A 

(T=t,4, cr=t,4- k,q 



With that, the electron-phonon interaction Hamiltonian (31) becomes 



7j Sq,7 4+q,aCk,<x iq,7 , (34) 



V-^ <7=t,; k q,7 
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where A^^^ = (b^^y + ) and the electron-phonon coupling constant 



The renormalized phonon frequencies can be found by determining the phonon propagator. 
We work in the random phase approximation (RPA), where the phonon propagator is given by 

f ■ \ r'i°^(q,iu;„) 
V^[q,iUn) = TTTs wr-- — r • (36) 

i i^T, {q.,iUJn) |gq| e(q,ia;„) 



In (36), 



2 (hiu^^l, 

(q, iu;.) = ^ \ , . (37) 



is the non-interacting phonon Matsubara Green function, 

1 1 

A /3 



n(q, iu;„) = WY.Y.^^ [^^°^(p + ^' + i^") ^^°^(P' (^8) 



iS7„ p 

is the the RPA electron polarization function and 

e(q, iun) = 1 - t^c(q ) n(q, iun) (39) 

is the RPA dielectric function. fc(q) = 2e^\ two-dimensional Fourier transform of the 
electron-electron Coulomb interaction potential fc(r) = ^^^^|^| . 



In (38) the non-interacting electron Green's function is defined as 



G^Z'^WiuJ-) = dre'-- {TrC^Ar)clA^))^ , (40) 

where (3 ^ ksT and Tr acts on the spin degrees of freedom a, a' =t, |. 



With ( 37 ), the RPA phonon propagator ( 36 ) can be written as 



VJci,iUn) = TTT, ^ TTT—^- (41) 

(ic^n)^-(^a^-2(/k.S)ig,/""'" 



7 1 £(q,iu;„) 
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The renormalized phonon frequency Wq,^ is obtained by seeking the poles of (41 1, which then 
gives 

i^o^nf = (K'li)^ + 2(H':i)lg.,.P . (42) 

where the analytic continuation iun to + iO+ is performed in order to obtain the retarded 



response functions. With (35 ), the RPA phonon propagator can finally be cast in the form 

h"^ ,on(q, 



{huj, 



(0)- 



+ 



where 21 is the primitive cell area. 



Ef 



Fermi 
Surface 



Intra-band 
Transition 




Inter-band 
Transition 



(43) 



Figure 14: Intra- and Inter-band transitions of DFQs that contribute to the renormalization of the 
prominent surface phonon branch, q is the phonon wave vector. 



The polarization function (38 ) can be conveniently decomposed into two contributions, one 
due to the intra-band and the other due to inter-band excitations as shown in figure |4] 



n(q, iUn) = nintra(q, l^n) + Ilinterlq, iUJn 



Upon performing the Matsubara sums in ( [38] ), we obtain 

'|k|<fcA 



n 



intra (.q? 



rp intra 

(2^ ^'"^ 



(C+q) - (C) ^ (^k+q) - ^F (^k ) 



k+q 
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k+q 



^k 



(44) 



(45) 



and 

ninter(q, l^n) 



|k|<fcA 



(27r 



|2 k,q 



(C+q) - (^k) , (^k+q) - (C) 



k+q 



k+q 



. (46) 



/ca denotes a momentum cutoff indicating the penetration of the metallic surface states into the 



bulk states, npiO 



e^^) ^ is the Fermi-Dirac distribution and, the chirality factors are 



77' intra 
-^k,q 



jp inter 
k,q 



1 / k-(k + q) 

2 V |k||k + q| 



(47) 



1 



1 



k • (k + q) 



2 V |k||k + q 
For T = and /i > the polarization expressions take the form 



n'"'-(q,ia;„ 

and 

n'°^-(q,i(^„) = 



|k|<fcA 



d^k 

(27r)2 



^(k, q) 



e(/i - Vf\^^ + q|) - 9(/x - VfM) 
Wirdk + ql - |k|) -iw„ 



(48) 



|k|<fcA 



d^k 

(2^ 



(k,q) 



e(;u-i;i.|k + q|) -1 ^ Q{^i-VF\k.\) - I 



fir(|k + q| + |k|) - iwn fF(|k + q| + |k| 



In order to solve (43 1 self consistently, we need to specify the momentum dependence of the 
coupling function Aq^. We argue that in the range of |q| < 2k p, the ionic screened potential 
y (q) has very weak dependence on g, since 2kp < gxF = 4^^^ = 0.5 A^^ the Thomas-Fermi 
wave vector. Moreover, invoking the sagittal-plane symmetry classification of the surface 



phonon modes, we specify the general polarization vector of the even parity mode as 

e^(q) = ei(q) + e||(q) , 

e±(q)~z, e||(q)~q. 



(49) 



The component (q) accounts for the vertical sheer component (ionic oscillations in a direc- 
tion perpendicular to the plane), while ey (q) is the longitudinal component. 



'The sagittal -plane is defined by the surface wave vector q and the surface normal. The only symmetry oper- 
ation associated with q along high-symmetry directions is a reflection through the sagittal-plane, having even and 
odd parity irreducible representations. 
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The corresponding coupling constant then reads 

Aq,7 = Aq • e^(q) = Aq • e_L(q) + Aq • e||(q) 

= Ax(q) + A||(q) . (50) 

In view of the near constancy of V{q) for q < 2kp, and the fact that the electron-phonon 
coupling involves the gradient of the screened potential, we write 

Aq = Af + Ma[°) (51) 

where go will be appropriately chosen as 2kF- The Dyson equation now becomes 

to first order in q. 

The model depends on two parameters: the bare phonon frequency and the Fermi wave 
vector. We identified the former with the experimental value of uj{q = 0) = 1.8 THz, where 
the DFQ response vanishes, kp was derived from a sample carrier concentration of —1.9 x 
10^^/cm'^, obtained from Hall measurements, which correspond to a Fermi energy of about 
300 meV and a Fermi wavevector of kp = 0.1 A^^, which is consistent with previous values 
reported for photoemission measurements. The coupling parameters X^f' and Ay^^ were left 
as fitting variables. Solutions for a;(q) were obtained by carrying out the integrals for 11 (q). 



and solving Eq. (43) iteratively. The calculated best-fit dispersion curve is shown in figure [T 



superimposed on the experimental data. 
The best fit parameters are 



(A^)2 = 10^ (meV)' ■ A , = 0.65 (53) 



M2t' ^ ' ' ' A 



with = 4x10 ^ meV, we obtain 



X^ = 50eV-A 
33 




this yields a real space value of 3.4 eV/A, which is quite reasonable, since it falls within the 
range of typical Coulombic interactions. 

We should note that recently, Thalmeier [|33l proposed a model based on acoustic -phonon 
perturbations, derived from the Dirac fermion Hamiltonian, but these perturbations are too weak 
to account for the observed dispersion. Moreover, as we demonstrated above, a TI actually 
presents a composite system consisting of an "ultra-thin metallic film" and an underlying in- 
sulating substrate. As such, the phonons in the metallic film can hardly penetrate into the 
substrate bulk in order to establish an acoustic Rayleigh branch. However, the unique feature 
of this system is that for q > 2kp, as DFQ screening becomes gradually suppressed, the sur- 
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face film becomes almost homogeneous with the underlying substrate, and the corresponding 
modes gradually penetrate into the bulk, ushering the V-shaped dispersion. We designate this 
unique behavior as a strong Kohn anomaly. Moreover, we should emphasize here that theo- 
retical modeUng and analysis of the interaction between the surface phonons and the DFQs 
demonstrate that the linear dispersion and isotropy of the DFQs are responsible for the profile 
of the prominent surface phonon branch for g < 2kF. 
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